R Notebook

library("rgdal") 
library("sf")
library('tidyverse')
library('spocc')
library('mapr')
library("scrubr")
library("osmdata")
library("raster")
library("osmdata")
library("taxize")
library("INLA")
library("data.table")
library("sp")
library("spdep")

Where and Why are people generating iNaturalist data in Boulder County?

iNaturalist is a mobile application that allows people to take photos of plants or animals they encounter in their environment and returns the identification of that species so people can learn more about those species. The user takes a photo of any species they encounter and upload that photo to the cloud, along with the metadata for the time and location when the photo was taken. Cloud computation services then use an AI classification algorithm to identify a list of candidate species matching the photo and location and return those result to the user for verification. The verified observation is then added to an iNaturalist databases that already includes over 68 Million individual occurrence records after only a few years of operation. These records are available through an iNaturalist API, but those records are also dumped into the GBIF database, where they undergo a small round of quality control that clean up the data in a few minor ways. I access these data using the ‘spocc’ package in R, which consolidates the API protocols for several different species occurrence API services into a single API syntax housed in one package. The ‘spocc’ package also interacts natively with the ‘taxis’ package to clean and standardize species names and common names across all taxa.

iNaturalist API

This code is an example showing how to pull data for a single species from the GBIF database. This is a great way to get data for an individual species, but GBIF suggests you use their online tool, rather than their API, for downloading multiple species at one time. You can set up a for loop to automate a list of species, but that unnecessarily drives up maintenance costs for GBIF, and they will again point you back to the online tool.

(df <- occ(query = 'Sciurus carolinensis', from = 'gbif',  has_coords = TRUE))
df2 <- occ2df(df)
class(df2)
df3 <- dframe(df2) %>% scrubr::fix_names( how = 'shortest')

df3 <- df3[df3$latitude > 45  & df3$longitude > -74 & df3$longitude < -73 ,,drop=TRUE]
df3

GBIF online download tool

Here I upload the file that I downloaded from gbif.org. To create this file, I used GBIF’s online tool to draw a bounding box around Boulder county and specified that I wanted all data with (1) Geographic coordinates, (2) No known errors, (3) recorded in 2021, and (4) iNaturalist listed as their provider. The provide that data in the form of a CSV, which I added to a folder called “rawData” in the working directory. These data come as a dataframe of occurrence records, which each record taking a row and different variables taking the columns. The first thing I do is clean the names using the ‘scrubr’ package. Then I plot the raw data to make sure they meet my expectations. I’m checking for anomalies that would be cause by improper loading, like the wrong spatial extent or unexpectedly empty fields.

gbif <- fread("rawData/0289344-200613084148143.csv",header=TRUE)
GBIF <- gbif %>% scrubr::fix_names( how = 'shortest')

plot(gbif$decimalLatitude, gbif$decimalLongitude)

I did a quick automated data clean using the ‘scrubr’ package, but I wanted to go a few steps further to replace scientific names with common names when appropriate. Common names are common in some fields like Ornithology and Mammalogy, but they are uncommon in Botany and Microbiology. The corrections produced in this section reflect those customs and only favor common names for species where they’re regularly provided in the iNaturalist dataset. I thinned the dataset to include only unique entries. This was largely about making this example smaller and easier to run quickly. There are repeated records from some notable individuals like birds of prey in public parks or notable trees on people’s walk to work. These individuals are recorded regularly and repeatedly and could be included in an example less pressed from space.

Hexagon grid to aggregate data for standardized analysis

Sampling scheme is a critical component to any analysis. I sample the points above by laying a hexagonal grid over the occurrence points and aggregating those points into those hex bins for modeling. For a 2D spatial analysis like the one we’re doing here, I prefer a hexagon lattice over a square lattice because it has a more symmetrical adjacency matrix. In a square lattice, the diagonals between cells are longer than the horizontal and vertical distances. In hexagonal grids, all adjacency distances are equal. This code creates the hex grid to be used below for sampling. The extent of this grid is set by the Boulder County bounding box and he resolution was set to be rather course to speed up computation time and keep this example light and fast.

# Define bounding box
 ext <- as(extent(getbb ("Boulder County Colorado")) , "SpatialPolygons")

# Set spatial projection
  crs(ext) <- "EPSG:4326"
  
# Use spsample to measure out a grid of center points within the bounding box  
  h <- spsample(ext, type = "hexagonal", cellsize = 0.01)
    
# convert center points to hexagons
  g <- HexPoints2SpatialPolygons(h, dx = 0.01)
  
# Reformat for easy plotting
  g <- st_as_sf(g)
 
# Plot to check that it's built corrrectly  
  plot(g, lwd=0.1)

  
# Make the dataset more presentable
  hex_polys <- cbind(seq(1, length(g$geometry)), g)
  colnames(hex_polys) <- c("id_polygons", "geometry") # change colnames
  

Aggregate iNaturalist data into hex bins

Now that we have iNaturalist data formatted as points and a hex grid built for sampling those points, it’s time to bin those points into their appropriate hexagons. I do this using an intersect function that builds a dataframe by assigning each occurrence point to the hexagonal polygon that encompasses it’s coordinates. I then group those data by their new polygon identifier and plot the hexagon grid with a fill color gradient representing the number of unique occurrence points in that polygon.

Inspect the points overlaying the grid.

st_crs(st_GBIF) <- crs(hex_polys)
intersection <- st_intersection(x = hex_polys, y = st_GBIF)

int_result <- intersection %>% 
  group_by(id_polygons) 

plot(g, lwd=0.1)
plot(int_result$geometry, pch=19, cex=0.25, col=adjustcolor("cornflowerblue", alpha.f = 0.2), add=TRUE)

Aggrigate points into hex bins

species_richness <- rep(0, max(hex_polys$id_polygons))
hex_counts <- cbind(hex_polys, species_richness)

int_count <- intersection %>% 
  group_by(id_polygons, .drop = FALSE) %>% 
  count()

hex_counts$species_richness[int_count$id_polygons] <- int_count$n

plot(hex_counts$geometry, lwd=0.001, 
     col=grey.colors(max(hex_counts$species_richness), rev = TRUE, start=0, end=1)[hex_counts$species_richness+1])
plot(int_result$geometry, pch=19, cex=0.05, col=adjustcolor("cornflowerblue", alpha.f = 0.2), add=TRUE)

Plot species richness using a fancy color scheme that works well for color-blind people.

A count of the number of species in an area is called species richness. Species richness is a common proxy variable used to describe differences between ecosystems or health within a single ecosystem. It provides an incomplete description of an ecosystem, but species richness is still frequently used as a good indicator of ecosystem type and health because more complex ecosystems have more species and so, the number of species can often approximate complexity and complexity is a defining feature of ecosystems.

plot(hex_counts$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness)+2, palette = "viridis", alpha = NULL, rev = TRUE, fixup = TRUE)[hex_counts$species_richness+1])

Adjacency matrix

We need to include a list of neighbor relationships as a random co-variate in our model. Her we calculate that adjacency matrix using a function that takes a list of polygons and returns the adjacency matrix. I give this function the hexagon grid, which is built of hexagonal polygons with centroids, and the function calculates the adjacency relationships for all those centroids.

hex_adj <- poly2nb(as(hex_counts, "Spatial") )

plot(hex_counts$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness)+2, palette = "viridis", alpha = NULL, rev = TRUE, fixup = TRUE)[hex_counts$species_richness+1])
plot(hex_adj, coordinates(as(hex_counts, "Spatial")), col="darkblue", lwd=0.2, add=TRUE)

Co-variates

Data don’t produce themselves in a vacuum, they are created by a process and that process involves many parts. I am now going to organize some data to act as co-variates in a model predicting iNaturalist species richness. These variables should relate to hypotheses about mechanisms that created these data and they should be at the same spatial and temporal scale as other data in the analysis.

Open Street Map data

On of my favorite datasets to use on a daily basis is Open Street Map. It is an amazing resource for a suite of different data about land use, amenities, administrative boundaries, points of interest, routing, and construction. The code below pulls polygons and multipolygons from the OSM API based on value/key pairs. I picked each value/key pair using the OSM tags wiki website, which catalogs all the available data in OSM. I search for different tags that described green space and I wrote a new API call for each value/key pair that I thought applicable. I concatenate those individual data requests into a single green space dataset. I then make three more individual calls to retrive data for buildings, roads, and agricultural fields. I think each of these may be contributing to our response variable so, I have left as individual dataframes so they can be treated individually in the modeling process.

name <- "Boulder County Colorado"

try(my_bbox <- osmdata::getbb (name))

my_bbox[1,1] <- my_bbox[1,1] + 0.02


dat1 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'park') %>%
    osmdata_sf ()
dat2 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'garden') %>%
    osmdata_sf ()
dat3 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'nature_reserve') %>%
    osmdata_sf ()
dat4 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'dog_park') %>%
    osmdata_sf ()
dat5 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'common') %>%
    osmdata_sf ()
dat6 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'pitch') %>%
    osmdata_sf ()
dat7 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'leisure', value = 'horse_riding') %>%
    osmdata_sf ()
dat8 <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'natural', value = 'wood') %>%
    osmdata_sf ()
greens <- c (dat1, dat2, dat3, dat4, dat5, dat6, dat7, dat8)


buildings <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'building') %>%
    osmdata_sf ()

roads <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'landuse', value = 'residential') %>%
    osmdata_sf ()

agriculture <- opq(bbox = my_bbox, timeout = 900) %>%
    add_osm_feature(key = 'landuse', value = 'farmland') %>%
    osmdata_sf()

Plot the green space data downloaded from OSM against the species richness data to make sure they overlap well.

plot(hex_counts$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness), palette = "viridis", alpha = 0, rev = TRUE, fixup = TRUE)[hex_counts$species_richness+1])
plot(greens$osm_polygons[1], add=TRUE, lwd=0.01, col=adjustcolor("grey", alpha.f = 1))

plot(greens$osm_multipolygons[1], add=TRUE, lwd=0.01, col=adjustcolor("grey", alpha.f = 1))
#plot(buildings$osm_polygons[1], add=TRUE, lwd=0.01, col=adjustcolor("black", alpha.f = 1))
plot(hex_counts$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness)+2, palette = "viridis", alpha = 0.8, rev = TRUE, fixup = TRUE)[hex_counts$species_richness+1], add=TRUE)

Aggregate covariates to hex bins using the same method as we did with species richness

Aggregating data to hexagon bins will be a required procedure for all covariates, which means we will need to repeat this procedure and it make sense to stop and make a function to keep our code clean. The first of these functions calculates polygon area for the greenspace polygons and adds that area calculation back into the dataset as a covariate.

covariate_hex_greens <- function(sf_object){

if(length(sf_object$osm_multipolygons) > 0){
  poly_area <- st_area(sf_object$osm_polygons)
multi_poly_area <- st_area(sf_object$osm_multipolygons)
log_area1 <- as.data.frame(as.numeric(log(poly_area)))
log_area2 <- as.data.frame(as.numeric(log(multi_poly_area)))
names(log_area1) <- "log_area"
names(log_area2) <- "log_area"
log_area <- rbind(log_area1, log_area2) 

coords <- rbind(sf_object$osm_polygons %>% st_centroid() %>% st_coordinates(), 
                sf_object$osm_multipolygons %>% st_centroid() %>% st_coordinates())
} else {
  poly_area <- st_area(sf_object$osm_polygons)
log_area <- as.data.frame(as.numeric(log(poly_area)))
names(log_area) <- "log_area"

coords <- sf_object$osm_polygons %>% st_centroid() %>% st_coordinates()
}

if(length(coords[,1]) != length(log_area[,1])) stop("wrong length")
area_points <- cbind(log_area, coords)
colnames(area_points)[2:3] <- c("longitude", "latitude")
sf_area_points <- st_as_sf(area_points, coords = c("longitude", "latitude"), 
    crs = 4326, agr = "constant")

return(sf_area_points)
}

The second function does the same aggregation procedure as I did with species richness, which is to say that it counts then number of occurrences within a bin and assigns that bin a value for that count. The other function in this chunk is a clone of the counting function modified to report mean rather than count. The mean version of the function is valuable or the greenspace area data where reporting the mean area of greenspace if more informative than the count of parks in that bin.

# Count the number of points in each bin
point_to_hex_count <- function(hex_polys, points){
intersection <- st_intersection(x = hex_polys, y = points)

counter <- rep(0, max(hex_polys$id_polygons))
hex_count <- cbind(hex_polys, counter)

counted <- intersection %>% 
  group_by(id_polygons, .drop = FALSE) %>% 
  count()

hex_count$area[counted$id_polygons] <- counted$n

return(hex_count)
}


# Report the mean value for each bin
point_to_hex_mean <- function(hex_polys, points){
intersection <- st_intersection(x = hex_polys, y = points)

area <- rep(0, max(hex_polys$id_polygons))
hex_area <- cbind(hex_polys, area)

avg_area <- intersection %>% 
  group_by(id_polygons, .drop = FALSE) %>% 
  summarise(mean = sum(log_area, na.rm=TRUE))

hex_area$area[avg_area$id_polygons] <- avg_area$mean

return(hex_area)
}

Run those functions on our green space data to produce a plot of green space area intensity

The plot below show the total area of greenspace per polygon. My hypothesis is that the presence of green space is a major predictor of iNaturalist data and want to compare this plot to the species richness plot to see it it looks like green spaces and iNaturalist records are generally located in the same places.

green_pts <- covariate_hex_greens(greens)
green_hex <- point_to_hex_mean(hex_polys, green_pts)
hex <- green_hex
plot_variable <- function(hex, variable ){
plot(hex$geometry, lwd=0.001, 
     col=hcl.colors(max(as.integer(green_hex$area)+2, na.rm = TRUE), palette = "viridis", alpha = NULL, rev = TRUE, fixup = TRUE)[as.integer(green_hex$area, na.rm=TRUE)+1])
}
plot_variable(green_hex, green_hex$area)

Combine Species richness and covariates into the same dataframe for analysis.

In the previous code, I have created two spatial objects that are both hex grids with aggregated data assigned to each bin. To combine those sets, I first drop the geometry from one of the two because the geometries are identical and they will show up as duplicates in the combined dataset. After droping the geometry from one, I do a left join so they match polygonID values and then everything shares the single remaining geometry column.

```r
hex_list <- st_drop_geometry(hex_counts) %>% inner_join(as.data.frame(hex_green_area), by=\id_polygons\) 
colnames(hex_list) <- c(\id_polygon\, \species_richness\, \park_size_average\,\geometry\)
hex_list

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



### INLA Bayesian inference


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuXG5oZXhfbGlzdCRsb2dfcGFya19hcmVhIDwtIGxvZyhoZXhfbGlzdCRwYXJrX3NpemVfYXZlcmFnZSlcbmhleF9saXN0W3doaWNoKGlzLmluZmluaXRlKGhleF9saXN0JGxvZ19wYXJrX2FyZWEpID09IFRSVUUpLFxcbG9nX3BhcmtfYXJlYVxcXSA8LSAwXG5cbiNSVzJkXG5tMC5ydzJkIDwtIGlubGEoc3BlY2llc19yaWNobmVzcyB+IGxvZ19wYXJrX2FyZWEgK1xuICAgIGYoaWRfcG9seWdvbiwgbW9kZWwgPSBcXHJ3MmRcXCwgbnJvdyA9NjQsIG5jb2wgPSA0MSksXG4gIGZhbWlseSA9IFxccG9pc3NvblxcLCBkYXRhID0gYXMuZGF0YS5mcmFtZShoZXhfbGlzdCksXG4gIGNvbnRyb2wucHJlZGljdG9yID0gbGlzdChjb21wdXRlID0gVFJVRSksXG4gIGNvbnRyb2wuY29tcHV0ZSA9IGxpc3QoZGljID0gVFJVRSkgKVxuXG5zdW1tYXJ5KG0wLnJ3MmQpXG5cbmhleF9saXN0JFJXMkQgPC0gbTAucncyZCRzdW1tYXJ5LmZpdHRlZC52YWx1ZXNbLCBcXDAuNXF1YW50XFxdXG5oZXhfbGlzdCRSVzJEc2QgPC0gbTAucncyZCRzdW1tYXJ5LmZpdHRlZC52YWx1ZXNbLCBcXHNkXFxdXG5gYGBcbmBgYCJ9 -->

```r
```r

hex_list$log_park_area <- log(hex_list$park_size_average)
hex_list[which(is.infinite(hex_list$log_park_area) == TRUE),\log_park_area\] <- 0

#RW2d
m0.rw2d <- inla(species_richness ~ log_park_area +
    f(id_polygon, model = \rw2d\, nrow =64, ncol = 41),
  family = \poisson\, data = as.data.frame(hex_list),
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE) )

summary(m0.rw2d)

hex_list$RW2D <- m0.rw2d$summary.fitted.values[, \0.5quant\]
hex_list$RW2Dsd <- m0.rw2d$summary.fitted.values[, \sd\]

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->




### Plot the mean predicted values from the model. 

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuI3Bsb3RfdmFyaWFibGUoaGV4X2xpc3QsIGZsb29yKGhleF9saXN0JFJXMkQpKVxuXG5tb2RlbF9wcmVkaWN0IDwtIGFzLmludGVnZXIoaGV4X2xpc3QkUlcyRClcbm1vZGVsX3ByZWRpY3Rbd2hpY2goaXMubmEobW9kZWxfcHJlZGljdCkgPT0gVFJVRSldIDwtIDAgXG4jbW9kZWxfcHJlZGljdFt3aGljaChtb2RlbF9wcmVkaWN0IDwgMCldIDwtIDAgXG4jbW9kZWxfcHJlZGljdFt3aGljaChtb2RlbF9wcmVkaWN0ID4gMjApXSA8LSAwIFxuXG5wbG90KGhleF9jb3VudHMkZ2VvbWV0cnksIGx3ZD0wLjAwMSwgXG4gICAgIGNvbD1oY2wuY29sb3JzKG1heChoZXhfY291bnRzJHNwZWNpZXNfcmljaG5lc3MpKzIsIHBhbGV0dGUgPSBcXHZpcmlkaXNcXCwgYWxwaGEgPSBOVUxMLCByZXYgPSBUUlVFLCBmaXh1cCA9IEZBTFNFKVtoZXhfY291bnRzJHNwZWNpZXNfcmljaG5lc3MrMV0pXG5cbnBsb3QoaGV4X2xpc3QkZ2VvbWV0cnksIGx3ZD0wLjAwMSwgXG4gICAgIGNvbD1oY2wuY29sb3JzKG1heChoZXhfY291bnRzJHNwZWNpZXNfcmljaG5lc3MpKzIsIHBhbGV0dGUgPSBcXHZpcmlkaXNcXCwgYWxwaGEgPSBOVUxMLCByZXYgPSBUUlVFLCBmaXh1cCA9IEZBTFNFKVttb2RlbF9wcmVkaWN0KzFdKVxuI3Bsb3QoaW50X3Jlc3VsdCRnZW9tZXRyeSwgcGNoPTE5LCBjZXg9MC4wNSwgY29sPWFkanVzdGNvbG9yKFxcY29ybmZsb3dlcmJsdWVcXCwgYWxwaGEuZiA9IDAuMiksIGFkZD1UUlVFKVxuYGBgXG5gYGAifQ== -->

```r
```r
#plot_variable(hex_list, floor(hex_list$RW2D))

model_predict <- as.integer(hex_list$RW2D)
model_predict[which(is.na(model_predict) == TRUE)] <- 0 
#model_predict[which(model_predict < 0)] <- 0 
#model_predict[which(model_predict > 20)] <- 0 

plot(hex_counts$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness)+2, palette = \viridis\, alpha = NULL, rev = TRUE, fixup = FALSE)[hex_counts$species_richness+1])

plot(hex_list$geometry, lwd=0.001, 
     col=hcl.colors(max(hex_counts$species_richness)+2, palette = \viridis\, alpha = NULL, rev = TRUE, fixup = FALSE)[model_predict+1])
#plot(int_result$geometry, pch=19, cex=0.05, col=adjustcolor(\cornflowerblue\, alpha.f = 0.2), add=TRUE)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->





<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYmVpLnRyZWVzMiRSVzJEIDwtIG0wLnJ3MmQkc3VtbWFyeS5maXR0ZWQudmFsdWVzWywgXFxtZWFuXFxdXG5cbmJlaS50cmVlczIkTUFURVJOMkQgPC0gbTAubTJkJHN1bW1hcnkuZml0dGVkLnZhbHVlc1ssIFxcbWVhblxcXVxuYGBgXG5gYGAifQ== -->

```r
```r
bei.trees2$RW2D <- m0.rw2d$summary.fitted.values[, \mean\]

bei.trees2$MATERN2D <- m0.m2d$summary.fitted.values[, \mean\]

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->




<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZ19zZiA8LSBzdF9hc19zZihnKVxuSFEgPC0gcmFzdGVyOjpleHRyYWN0KHJhc3Rlcl90aGlzX0hRLCAgZ19zZiAsIGZ1biA9IHN1bSwgbmEucm0gPSBUUlVFLCBkZiA9IFRSVUUsIG5vcm1hbGl6ZVdlaWdodHM9VFJVRSkgIFxuYGBgXG5gYGAifQ== -->

```r
```r
g_sf <- st_as_sf(g)
HQ <- raster::extract(raster_this_HQ,  g_sf , fun = sum, na.rm = TRUE, df = TRUE, normalizeWeights=TRUE)  

```